##########################################################
### Disorganized Political Violence:                   ###
### A Demonstration Case of Temperature and Insurgency ###
### Replication Code: Primary Results --               ###
### Temperature Response Function Results &            ###
### Comparisons (Appendix Results)                     ###
### Andrew Shaver, Alex Bollfrass                      ###
### Date: 07/04/2022                                   ###
##########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(arm)
library(broom)
library(dplyr)
library(ggpubr)
library(ggplot2)
library(lfe)
library(lubridate)
library(MASS)
library(mgcv)
library(stargazer)
library(sandwich)
library(nnet)
library(scales)
library(gridExtra)
library(grid)

packageVersion("arm")
# ‘1.12.2’
packageVersion("broom")
# ‘0.7.12’
packageVersion("dplyr")
# ‘1.0.9’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("ggplot2")
# ‘3.3.6’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("lubridate")
# ‘1.8.0’
packageVersion("MASS")
# ‘7.3.55’
packageVersion("mgcv")
# ‘1.8.39’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("sandwich")
# ‘3.0.1’
packageVersion("nnet")
# ‘7.3.17’
packageVersion("scales")
# ‘1.2.0’
packageVersion("gridExtra")
# ‘2.3’
packageVersion("grid")
# ‘4.1.2’

# Load data:
# For primary temperature-violence analysis:
# Daily analysis:
iq_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iraq_baghdad_basrah_processed_data.rds")
af_sig1 <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/afghanstan_viol_districts_processed_data.rds")
# For attitude analysis:
bg_sur <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/baghdad_survey.rds")
# For median analysis:
iq_median_temp <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/Affective_Motivations/MajorRevision/R3_Submission/ProcessedData/iq_daily_median.rds")

# Construct range of all temperatures observed across Afghanistan and Iraq:
all_observed_temps <- c(iq_sig1$max, af_sig1$max)
all_observed_temps <- as.data.frame(all_observed_temps)
colnames(all_observed_temps) <- "temperature"

# For temperature deviations for the Baghdad survey data, given that the unit is not the
# match deviations from panel:
# Columns to keep:
to_keep <- c("date", "temp_1", "temp_2", "temp_3", "temp_4", 
             "temp_5", "temp_6", "temp_7", "daylight")
bg_sur <- plyr::join(bg_sur, iq_sig1[iq_sig1$governorate %in% "Baghdad", to_keep], by = "date")

# Create survey dataset removing top income earners:
bg_sur.s3.l <- bg_sur[bg_sur$income.w < 75000, ]

# Find minimum and maximum temperatures:
temperature_range <- round(min(iq_sig1$max, af_sig1$max, bg_sur$max)):round(max(iq_sig1$max, af_sig1$max, bg_sur$max))

###
### IRAQ.
###

# Baghdad and Basrah:
# Unconstrained Violence:
# OLS model results:

IQ.w.u.l <- felm(log(unconstrained+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.q <- felm(log(unconstrained+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.u.c <- felm(log(unconstrained+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)

IQ.m.u.l <- felm(log(unconstrained+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.u.q <- felm(log(unconstrained+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.u.c <- felm(log(unconstrained+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)

# Constrained Violence:
IQ.w.c.l <- felm(log(constrained+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.c.q <- felm(log(constrained+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.c.c <- felm(log(constrained+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)

IQ.m.c.l <- felm(log(constrained+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.c.q <- felm(log(constrained+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.c.c <- felm(log(constrained+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)

# IED:
IQ.w.i.l <- felm(log(ied+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.i.q <- felm(log(ied+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)
IQ.w.i.c <- felm(log(ied+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + governorate | 0 | 0, data = iq_sig1)

IQ.m.i.l <- felm(log(ied+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.i.q <- felm(log(ied+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)
IQ.m.i.c <- felm(log(ied+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                   unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + governorate | 0 | 0, data = iq_sig1)

###
### AFGHANISTAN.
###

AF.w.u.l <- felm(log(df+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.u.q <- felm(log(df+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.u.c <- felm(log(df+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)

AF.m.u.l <- felm(log(df+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.u.q <- felm(log(df+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.u.c <- felm(log(df+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)

AF.w.i.l <- felm(log(ied+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.i.q <- felm(log(ied+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)
AF.w.i.c <- felm(log(ied+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   week + station | 0 | month, data = af_sig1)

AF.m.i.l <- felm(log(ied+1) ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.i.q <- felm(log(ied+1) ~ max + max2 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)
AF.m.i.c <- felm(log(ied+1) ~ max + max2 + max3 + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                   df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                   idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                   ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                   month + station | 0 | month, data = af_sig1)

###
### Create dataset for prediction:
###

new_data_iq <- data.frame(max = temperature_range, 
                          max2 = (temperature_range)^2,
                          max3 = (temperature_range)^3,
                          
                          prcp = mean(iq_sig1$prcp, na.rm = T), 
                          wdsp = mean(iq_sig1$wdsp, na.rm = T), 
                          dewp = mean(iq_sig1$dewp, na.rm = T), 
                          visib = mean(iq_sig1$visib, na.rm = T), 
                          mxspd = mean(iq_sig1$mxspd, na.rm = T), 
                          daylight = mean(iq_sig1$daylight, na.rm = T), 
                          hoursofpower = mean(iq_sig1$hoursofpower, na.rm = T), 
                          
                          unconstrained = mean(iq_sig1$unconstrained), 
                          constrained = mean(iq_sig1$constrained),  
                          idf = mean(iq_sig1$idf),  
                          ied = mean(iq_sig1$ied),  
                          
                          temp_1 = mean(iq_sig1$temp_1, na.rm = T),
                          temp_2 = mean(iq_sig1$temp_2, na.rm = T),
                          temp_3 = mean(iq_sig1$temp_3, na.rm = T),
                          temp_4 = mean(iq_sig1$temp_4, na.rm = T),
                          temp_5 = mean(iq_sig1$temp_5, na.rm = T),
                          temp_6 = mean(iq_sig1$temp_6, na.rm = T),
                          temp_7 = mean(iq_sig1$temp_7, na.rm = T),
                          
                          unconstrained_1 = mean(iq_sig1$unconstrained_1, na.rm = T),
                          unconstrained_2 = mean(iq_sig1$unconstrained_2, na.rm = T),
                          unconstrained_3 = mean(iq_sig1$unconstrained_3, na.rm = T),
                          unconstrained_4 = mean(iq_sig1$unconstrained_4, na.rm = T),
                          unconstrained_5 = mean(iq_sig1$unconstrained_5, na.rm = T),
                          unconstrained_6 = mean(iq_sig1$unconstrained_6, na.rm = T),
                          unconstrained_7 = mean(iq_sig1$unconstrained_7, na.rm = T),
                          
                          constrained_1 = mean(iq_sig1$constrained_1, na.rm = T),
                          constrained_2 = mean(iq_sig1$constrained_2, na.rm = T),
                          constrained_3 = mean(iq_sig1$constrained_3, na.rm = T),
                          constrained_4 = mean(iq_sig1$constrained_4, na.rm = T),
                          constrained_5 = mean(iq_sig1$constrained_5, na.rm = T),
                          constrained_6 = mean(iq_sig1$constrained_6, na.rm = T),
                          constrained_7 = mean(iq_sig1$constrained_7, na.rm = T),
                          
                          idf_1 = mean(iq_sig1$idf_1, na.rm = T),
                          idf_2 = mean(iq_sig1$idf_2, na.rm = T),
                          idf_3 = mean(iq_sig1$idf_3, na.rm = T),
                          idf_4 = mean(iq_sig1$idf_4, na.rm = T),
                          idf_5 = mean(iq_sig1$idf_5, na.rm = T),
                          idf_6 = mean(iq_sig1$idf_6, na.rm = T),
                          idf_7 = mean(iq_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(iq_sig1$ied_1, na.rm = T),
                          ied_2 = mean(iq_sig1$ied_2, na.rm = T),
                          ied_3 = mean(iq_sig1$ied_3, na.rm = T),
                          ied_4 = mean(iq_sig1$ied_4, na.rm = T),
                          ied_5 = mean(iq_sig1$ied_5, na.rm = T),
                          ied_6 = mean(iq_sig1$ied_6, na.rm = T),
                          ied_7 = mean(iq_sig1$ied_7, na.rm = T),
                          
                          week = median(iq_sig1$week),
                          month = median(iq_sig1$month),
                          governorate = "Baghdad")

new_data_af <- data.frame(max = temperature_range, 
                          max2 = (temperature_range)^2,
                          max3 = (temperature_range)^3,
                          
                          prcp = mean(af_sig1$prcp, na.rm = T), 
                          wdsp = mean(af_sig1$wdsp, na.rm = T), 
                          dewp = mean(af_sig1$dewp, na.rm = T), 
                          visib = mean(af_sig1$visib, na.rm = T), 
                          mxspd = mean(af_sig1$mxspd, na.rm = T), 
                          daylight = mean(af_sig1$daylight, na.rm = T), 
                          
                          df = mean(af_sig1$df), 
                          idf = mean(af_sig1$idf),  
                          ied = mean(af_sig1$ied), 
                          
                          temp_1 = mean(af_sig1$temp_1, na.rm = T),
                          temp_2 = mean(af_sig1$temp_2, na.rm = T),
                          temp_3 = mean(af_sig1$temp_3, na.rm = T),
                          temp_4 = mean(af_sig1$temp_4, na.rm = T),
                          temp_5 = mean(af_sig1$temp_5, na.rm = T),
                          temp_6 = mean(af_sig1$temp_6, na.rm = T),
                          temp_7 = mean(af_sig1$temp_7, na.rm = T),
                          
                          df_1 = mean(af_sig1$df_1, na.rm = T),
                          df_2 = mean(af_sig1$df_2, na.rm = T),
                          df_3 = mean(af_sig1$df_3, na.rm = T),
                          df_4 = mean(af_sig1$df_4, na.rm = T),
                          df_5 = mean(af_sig1$df_5, na.rm = T),
                          df_6 = mean(af_sig1$df_6, na.rm = T),
                          df_7 = mean(af_sig1$df_7, na.rm = T),
                          
                          idf_1 = mean(af_sig1$idf_1, na.rm = T),
                          idf_2 = mean(af_sig1$idf_2, na.rm = T),
                          idf_3 = mean(af_sig1$idf_3, na.rm = T),
                          idf_4 = mean(af_sig1$idf_4, na.rm = T),
                          idf_5 = mean(af_sig1$idf_5, na.rm = T),
                          idf_6 = mean(af_sig1$idf_6, na.rm = T),
                          idf_7 = mean(af_sig1$idf_7, na.rm = T),
                          
                          ied_1 = mean(af_sig1$ied_1, na.rm = T),
                          ied_2 = mean(af_sig1$ied_2, na.rm = T),
                          ied_3 = mean(af_sig1$ied_3, na.rm = T),
                          ied_4 = mean(af_sig1$ied_4, na.rm = T),
                          ied_5 = mean(af_sig1$ied_5, na.rm = T),
                          ied_6 = mean(af_sig1$ied_6, na.rm = T),
                          ied_7 = mean(af_sig1$ied_7, na.rm = T),
                          
                          week = median(af_sig1$week),
                          month = median(af_sig1$month),
                          station = 409900)
###
### Then, get predictions:
###

# The following code is sourced from Stack Overflow:
# https://stackoverflow.com/questions/30491545/predict-method-for-felm-from-lfe-package

predict.felm <- function(object, newdata, se.fit = FALSE,
                         interval = "none",
                         level = 0.95){
  if(missing(newdata)){
    stop("predict.felm requires newdata and predicts for all group effects = 0.")
  }
  
  tt <- terms(object)
  Terms <- delete.response(tt)
  attr(Terms, "intercept") <- 0
  
  m.mat <- model.matrix(Terms, data = newdata)
  m.coef <- as.numeric(object$coef)
  fit <- as.vector(m.mat %*% object$coef)
  fit <- data.frame(fit = fit)
  
  if(se.fit | interval != "none"){
    if(!is.null(object$clustervcv)){
      vcov_mat <- object$clustervcv
    } else if (!is.null(object$robustvcv)) {
      vcov_mat <- object$robustvcv
    } else if (!is.null(object$vcv)){
      vcov_mat <- object$vcv
    } else {
      stop("No vcv attached to felm object.")
    }
    se.fit_mat <- sqrt(diag(m.mat %*% vcov_mat %*% t(m.mat)))
  }
  if(interval == "confidence"){
    t_val <- qt((1 - level) / 2 + level, df = object$df.residual)
    fit$lwr <- fit$fit - t_val * se.fit_mat
    fit$upr <- fit$fit + t_val * se.fit_mat
  } else if (interval == "prediction"){
    stop("interval = \"prediction\" not yet implemented")
  }
  if(se.fit){
    return(list(fit=fit, se.fit=se.fit_mat))
  } else {
    return(fit)
  }
}

pred_IQ.w.u.l <- predict(IQ.w.u.l, new_data_iq)
pred_IQ.w.u.q <- predict(IQ.w.u.q, new_data_iq)
pred_IQ.w.u.c <- predict(IQ.w.u.c, new_data_iq)

pred_IQ.m.u.l <- predict(IQ.m.u.l, new_data_iq)
pred_IQ.m.u.q <- predict(IQ.m.u.q, new_data_iq)
pred_IQ.m.u.c <- predict(IQ.m.u.c, new_data_iq)

pred_IQ.w.c.l <- predict(IQ.w.c.l, new_data_iq)
pred_IQ.w.c.q <- predict(IQ.w.c.q, new_data_iq)
pred_IQ.w.c.c <- predict(IQ.w.c.c, new_data_iq)

pred_IQ.m.c.l <- predict(IQ.m.c.l, new_data_iq)
pred_IQ.m.c.q <- predict(IQ.m.c.q, new_data_iq)
pred_IQ.m.c.c <- predict(IQ.m.c.c, new_data_iq)

pred_IQ.w.i.l <- predict(IQ.w.i.l, new_data_iq)
pred_IQ.w.i.q <- predict(IQ.w.i.q, new_data_iq)
pred_IQ.w.i.c <- predict(IQ.w.i.c, new_data_iq)

pred_IQ.m.i.l <- predict(IQ.m.i.l, new_data_iq)
pred_IQ.m.i.q <- predict(IQ.m.i.q, new_data_iq)
pred_IQ.m.i.c <- predict(IQ.m.i.c, new_data_iq)

# Afghanistan:
pred_AF.w.u.l <- predict(AF.w.u.l, new_data_af)
pred_AF.w.u.q <- predict(AF.w.u.q, new_data_af)
pred_AF.w.u.c <- predict(AF.w.u.c, new_data_af)

pred_AF.m.u.l <- predict(AF.m.u.l, new_data_af)
pred_AF.m.u.q <- predict(AF.m.u.q, new_data_af)
pred_AF.m.u.c <- predict(AF.m.u.c, new_data_af)

pred_AF.w.i.l <- predict(AF.w.i.l, new_data_af)
pred_AF.w.i.q <- predict(AF.w.i.q, new_data_af)
pred_AF.w.i.c <- predict(AF.w.i.c, new_data_af)

pred_AF.m.i.l <- predict(AF.m.i.l, new_data_af)
pred_AF.m.i.q <- predict(AF.m.i.q, new_data_af)
pred_AF.m.i.c <- predict(AF.m.i.c, new_data_af)

###
### Then, build dataset from the predictions (for statistically significant models only:
###

preds_ols <- cbind(pred_IQ.w.u.q, 
                   pred_IQ.w.u.c,
                   pred_IQ.m.u.l,
                   pred_IQ.m.u.q, 
                   pred_AF.w.u.l, pred_AF.w.u.q, pred_AF.w.u.c, 
                   pred_AF.m.u.q, pred_AF.m.u.c)
colnames(preds_ols) <- c("model1_iq1", "model2_iq2", 
                         "model3_iq3", "model4_iq4",
                      
                         "model5_af1", "model6_af2", "model7_af3", 
                         "model8_af4", "model9_af5")
# Create temperature range:
preds_ols$temperature <- temperature_range

# Take the means:
preds_ols$mean_iq <- rowMeans(preds_ols[, 1:2])
preds_ols$mean_af <- rowMeans(preds_ols[, 3:8])

# Then, distinguish predictions for temperatures outside of each country's observed values:
temperature_range_iq <- round(min(iq_sig1$max)):round(max(iq_sig1$max))
temperature_range_af <- round(min(af_sig1$max)):round(max(af_sig1$max))

preds_ols$model1_iq1_c <- 
  preds_ols$model2_iq2_c <- 
  preds_ols$model3_iq3_c <- 
  preds_ols$model4_iq4_c <- 
  preds_ols$model5_af1_c <- preds_ols$model6_af2_c <- preds_ols$model7_af3_c <- 
  preds_ols$model8_af4_c <- preds_ols$model9_af5_c <- NA

preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model1_iq1_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model1_iq1
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model2_iq2_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model2_iq2
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model3_iq3_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model3_iq3
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model4_iq4_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model4_iq4

preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model5_af1_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model5_af1
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model6_af2_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model6_af2
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model7_af3_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model7_af3
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model8_af4_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model8_af4
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model9_af5_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model9_af5

preds_ols$mean_iq_c <- NA
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$mean_iq_c <- as.vector(preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$mean_iq)
preds_ols$mean_af_c <- NA
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$mean_af_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$mean_af

###
### Then, plot:
###

p_combined <- ggplot() + 

  geom_line(data = preds_ols, aes(x = temperature, y = model1_iq1), color = "orange3", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model2_iq2), color = "orange3", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model3_iq3), color = "orange3", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model4_iq4), color = "orange3", size = 0.15, linetype = "dotted") +
  
  geom_line(data = preds_ols, aes(x = temperature, y = model5_af1), color = "deepskyblue1", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model6_af2), color = "deepskyblue1", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model7_af3), color = "deepskyblue1", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model8_af4), color = "deepskyblue1", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model9_af5), color = "deepskyblue1", size = 0.15, linetype = "dashed") +
  
  geom_line(data = preds_ols, aes(x = temperature, y = model1_iq1_c), color = "orange3", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model2_iq2_c), color = "orange3", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model3_iq3_c), color = "orange3", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model4_iq4_c), color = "orange3", size = 0.5, linetype = "dotted") +
  
  geom_line(data = preds_ols, aes(x = temperature, y = model5_af1_c), color = "deepskyblue1", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model6_af2_c), color = "deepskyblue1", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model7_af3_c), color = "deepskyblue1", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model8_af4_c), color = "deepskyblue1", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model9_af5_c), color = "deepskyblue1", size = 0.5, linetype = "dashed") +
  
  geom_line(data = preds_ols, aes(x = temperature, y = mean_iq_c), color = "orange3", size = 1) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_af_c), color = "deepskyblue1", size = 1) +
  
  geom_line(data = preds_ols, aes(x = temperature, y = mean_iq), color = "orange3", size = 0.5) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_af), color = "deepskyblue1", size = 0.5) +
  
  theme_bw() +
  geom_vline(xintercept = preds_ols[which(preds_ols$mean_iq == max(preds_ols$mean_iq)), ]$temperature, linetype=4, color = "orange3") +
  geom_vline(xintercept = preds_ols[which(preds_ols$mean_af == max(preds_ols$mean_af)), ]$temperature-0.25, linetype=4, color = "deepskyblue1") +
  
  ggtitle("") +
  ylab("Predicted Instances of Least Constrained Violent Attacks Across Models") + 
  xlab("Temperature (°F)") +
  xlim(0,122) + 
  ylim(0,5.25) + 
  
  scale_x_continuous(breaks = seq(9, 122, by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=18),
    axis.text.x = element_text(angle = 45)) + 
  annotate("segment", x = 9, xend = 15, y = 5, yend = 5, colour = "deepskyblue1", size=1, alpha=0.6) +
  annotate(geom="text", x=16, y= 5, label="Afghanistan",
           color="black", size = 4, , hjust = 0) +
  annotate("segment", x = 9, xend = 15, y = 4.5, yend = 4.5, colour = "orange3", size=1, alpha=0.6) +
  annotate(geom="text", x=16, y= 4.5, label="Iraq",
           color="black", size = 4, , hjust = 0) +
  annotate("segment", x = 9, xend = 15, y = 4.5, yend = 4.5, colour = "orange3", size=1, alpha=0.6) +
  
  annotate("segment", x = 31.5, xend = 37.5, y = 5, yend = 5, colour = "black", 
           size=0.75, alpha=0.6, linetype = "dashed") +
  annotate("segment", x = 31.5, xend = 37.5, y = 4.5, yend = 4.5, colour = "black",
           size=0.75, alpha=0.6, linetype = "dotted") +
  
  annotate(geom="text", x=38, y= 5, label="Week FE models",
           color="black", size = 4, , hjust = 0) +
  annotate(geom="text", x=38, y= 4.5, label="Month FE models",
           color="black", size = 4, , hjust = 0) +
  
  geom_rect(aes(xmin = 8, xmax = 60, ymin = 4.25, ymax = 5.25),
            fill = "transparent", color = "black", size = 0.5) +
  
  annotate("segment", x = 40, xend = preds_ols[which(preds_ols$mean_iq == max(preds_ols$mean_iq)), ]$temperature, 
           y = 4, yend = max(preds_ols$mean_iq), colour = "gray", size=0.75, alpha=0.6) +
  annotate("segment", x = 40, xend = preds_ols[which(preds_ols$mean_af == max(preds_ols$mean_af)), ]$temperature, 
           y = 4, yend = max(preds_ols$mean_af), colour = "gray", size=0.75, alpha=0.6) +
  annotate(geom="text", x=39, y= 4, label="°F at which violence at max",
           color="black", size = 4, , hjust = 1)

p_combined

# Next, generate predicted counts:

###
### IRAQ.
###

# Baghdad and Basrah:
# Unconstrained Violence:
# QP model results:
IQ.w.u.l.qp <- bayesglm(unconstrained ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.u.q.qp <- bayesglm(unconstrained ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.u.c.qp <- bayesglm(unconstrained ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.u.l.qp <- bayesglm(unconstrained ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.u.q.qp <- bayesglm(unconstrained ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.u.c.qp <- bayesglm(unconstrained ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + constrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# Constrained Violence:
IQ.w.c.l.qp <- bayesglm(constrained ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.c.q.qp <- bayesglm(constrained ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.c.c.qp <- bayesglm(constrained ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.c.l.qp <- bayesglm(constrained ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.c.q.qp <- bayesglm(constrained ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.c.c.qp <- bayesglm(constrained ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + ied + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

# IED:
IQ.w.i.l.qp <- bayesglm(ied ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.i.q.qp <- bayesglm(ied ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.w.i.c.qp <- bayesglm(ied ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

IQ.m.i.l.qp <- bayesglm(ied ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.i.q.qp <- bayesglm(ied ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)
IQ.m.i.c.qp <- bayesglm(ied ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower + idf + constrained + unconstrained +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                          unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(governorate), family = "quasipoisson", data = iq_sig1)

###
### AFGHANISTAN.
###

AF.w.u.l.qp <- bayesglm(df ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.u.q.qp <- bayesglm(df ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.u.c.qp <- bayesglm(df ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)

AF.m.u.l.qp <- bayesglm(df ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.u.q.qp <- bayesglm(df ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.u.c.qp <- bayesglm(df ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + idf + ied +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)

AF.w.i.l.qp <- bayesglm(ied ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.i.q.qp <- bayesglm(ied ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.w.i.c.qp <- bayesglm(ied ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(week) + as.factor(district), family = "quasipoisson", data = af_sig1)

AF.m.i.l.qp <- bayesglm(ied ~ max + prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.i.q.qp <- bayesglm(ied ~ max + max2 +  prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)
AF.m.i.c.qp <- bayesglm(ied ~ max + max2 +  max3 +  prcp + wdsp + dewp + visib + mxspd + daylight + df + idf +
                          temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                          df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                          idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                          ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                          as.factor(month) + as.factor(district), family = "quasipoisson", data = af_sig1)

# Next, generate predicted counts:

# First, construct a function for generating counts and creating output that can be plotted.
plot_function <- function(reg_output, X, dataset, type, model){
  
  dataset <- dataset[!(is.na(dataset$max)), ]
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(reg_output), Sigma = vcov(reg_output))
  paste <- model.matrix(reg_output)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps.df1 <- matrix(NA, length(unique(round(dataset$max))), 100)
  for(i in 1:length(unique(round(dataset$max)))){
    if(X %in% "linear"){
      paste[, 2] <- i + (min(round(dataset$max), na.rm = T)-1)
      x.beta <- paste %*% t(beta.sim)
      if(model %in% "count"){
        yhat <- exp(x.beta)
      } else if(model %in% "logit"){
        yhat <- exp(x.beta) / (exp(x.beta) + 1)
      }
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "quadratic"){
      
      paste[, 2] <- i + (min(round(dataset$max), na.rm = T)-1)
      paste[, 3] <- (i + (min(round(dataset$max), na.rm = T)-1))^2
      x.beta <- paste %*% t(beta.sim)
      if(model %in% "count"){
        yhat <- exp(x.beta)
      } else if(model %in% "logit"){
        yhat <- exp(x.beta) / (exp(x.beta) + 1)
      }

      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(X %in% "cubic"){
      paste[, 2] <- i + (min(round(dataset$max), na.rm = T)-1)
      paste[, 3] <- (i + (min(round(dataset$max), na.rm = T)-1))^2
      paste[, 4] <- (i + (min(round(dataset$max), na.rm = T)-1))^3
      x.beta <- paste %*% t(beta.sim)
      if(model %in% "count"){
        yhat <- exp(x.beta)
      } else if(model %in% "logit"){
        yhat <- exp(x.beta) / (exp(x.beta) + 1)
      }
      reps.df1[i, ] <- apply(yhat, 2, mean, na.rm=T)
    }
  }
  
  plot_dataset <- as.data.frame(matrix(NA, length(unique(round(dataset$max))),3))
  colnames(plot_dataset) <- c("means", "cil", "ciu")
  
  plot_dataset$means <- apply(reps.df1, 1 , mean)
  se <- apply(reps.df1, 1 , sd)
  plot_dataset$cil <- plot_dataset$means + qnorm(0.975)*se
  plot_dataset$ciu <- plot_dataset$means - qnorm(0.975)*se
  plot_dataset$temperature <- sort(unique(round(dataset$max)))
  
  if(X %in% "linear"){
    plot_dataset$model <- "linear"
  } else if(X %in% "quadratic"){
    plot_dataset$model <- "quadratic"
  } else if(X %in% "cubic"){
    plot_dataset$model <- "cubic"
  }
  
  plot_dataset$max <- 0
  plot_dataset[which.max(plot_dataset$means), ]$max <- 1
  plot_dataset$attack_type <- type
  
  return(plot_dataset)
}

# Iraq:

p.IQ.w.u.l.qp <- plot_function(IQ.w.u.l.qp, "linear", iq_sig1, "unconstrained", "count")
p.IQ.w.u.q.qp <- plot_function(IQ.w.u.q.qp, "quadratic", iq_sig1, "unconstrained", "count")
p.IQ.w.u.c.qp <- plot_function(IQ.w.u.c.qp, "cubic", iq_sig1, "unconstrained", "count")

p.IQ.m.u.l.qp <- plot_function(IQ.m.u.l.qp, "linear", iq_sig1, "unconstrained", "count")
p.IQ.m.u.q.qp <- plot_function(IQ.m.u.q.qp, "quadratic", iq_sig1, "unconstrained", "count")
p.IQ.m.u.c.qp <- plot_function(IQ.m.u.c.qp, "cubic", iq_sig1, "unconstrained", "count")

p.IQ.w.c.l.qp <- plot_function(IQ.w.c.l.qp, "linear", iq_sig1, "constrained", "count")
p.IQ.w.c.q.qp <- plot_function(IQ.w.c.q.qp, "quadratic", iq_sig1, "constrained", "count")
p.IQ.w.c.c.qp <- plot_function(IQ.w.c.c.qp, "cubic", iq_sig1, "constrained", "count")

p.IQ.m.c.l.qp <- plot_function(IQ.m.c.l.qp, "linear", iq_sig1, "constrained", "count")
p.IQ.m.c.q.qp <- plot_function(IQ.m.c.q.qp, "quadratic", iq_sig1, "constrained", "count")
p.IQ.m.c.c.qp <- plot_function(IQ.m.c.c.qp, "cubic", iq_sig1, "constrained", "count")

p.IQ.w.i.l.qp <- plot_function(IQ.w.i.l.qp, "linear", iq_sig1, "ied", "count")
p.IQ.w.i.q.qp <- plot_function(IQ.w.i.q.qp, "quadratic", iq_sig1, "ied", "count")
p.IQ.w.i.c.qp <- plot_function(IQ.w.i.c.qp, "cubic", iq_sig1, "ied", "count")

p.IQ.m.i.l.qp <- plot_function(IQ.m.i.l.qp, "linear", iq_sig1, "ied", "count")
p.IQ.m.i.q.qp <- plot_function(IQ.m.i.q.qp, "quadratic", iq_sig1, "ied", "count")
p.IQ.m.i.c.qp <- plot_function(IQ.m.i.c.qp, "cubic", iq_sig1, "ied", "count")

# Afghanistan:

p.AF.w.u.l.qp <- plot_function(AF.w.u.l.qp, "linear", af_sig1, "unconstrained", "count")
p.AF.w.u.q.qp <- plot_function(AF.w.u.q.qp, "quadratic", af_sig1, "unconstrained", "count")
p.AF.w.u.c.qp <- plot_function(AF.w.u.c.qp, "cubic", af_sig1, "unconstrained", "count")

p.AF.m.u.l.qp <- plot_function(AF.m.u.l.qp, "linear", af_sig1, "unconstrained", "count")
p.AF.m.u.q.qp <- plot_function(AF.m.u.q.qp, "quadratic", af_sig1, "unconstrained", "count")
p.AF.m.u.c.qp <- plot_function(AF.m.u.c.qp, "cubic", af_sig1, "unconstrained", "count")

p.AF.w.i.l.qp <- plot_function(AF.w.i.l.qp, "linear", af_sig1, "ied", "count")
p.AF.w.i.q.qp <- plot_function(AF.w.i.q.qp, "quadratic", af_sig1, "ied", "count")
p.AF.w.i.c.qp <- plot_function(AF.w.i.c.qp, "cubic", af_sig1, "ied", "count")

p.AF.m.i.l.qp <- plot_function(AF.m.i.l.qp, "linear", af_sig1, "ied", "count")
p.AF.m.i.q.qp <- plot_function(AF.m.i.q.qp, "quadratic", af_sig1, "ied", "count")
p.AF.m.i.c.qp <- plot_function(AF.m.i.c.qp, "cubic", af_sig1, "ied", "count")

# Plot:

# Iraq:
# Least Constrained:
plot.IQ.m.u.l.qp <- ggplot() + 

  geom_line(data = p.IQ.w.u.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.w.u.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  
  geom_line(data = p.IQ.m.u.l.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.m.u.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.IQ.m.u.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_sig1$max)), max(round(iq_sig1$max)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(2, 9) +
  ylab("") + 
  xlab("") +
  ggtitle("Least Constrained")

# Most Constrained:
plot.IQ.m.c.l.qp <- ggplot() + 
  geom_line(data = p.IQ.w.c.l.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.IQ.w.c.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.w.c.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.c.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.w.c.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.c.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  
  geom_line(data = p.IQ.m.c.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.c.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  geom_line(data = p.IQ.m.c.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.c.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "green4") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_sig1$max)), max(round(iq_sig1$max)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 7.5) +
  ylab("") + 
  xlab("") +
  ggtitle("Most Constrained")

# IED:
plot.IQ.m.i.l.qp <- ggplot() + 

  geom_line(data = p.IQ.w.i.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.IQ.w.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  geom_line(data = p.IQ.w.i.c.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.IQ.w.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  
  geom_line(data = p.IQ.m.i.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.IQ.m.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  geom_line(data = p.IQ.m.i.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.IQ.m.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +

  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(iq_sig1$max)), max(round(iq_sig1$max)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(4, 13) +
  ylab("") + 
  xlab("") +
  ggtitle("IED Attacks")

# Afghanistan:
# DF:
plot.AF.m.u.l.qp <- ggplot() + 
  geom_line(data = p.AF.w.u.l.qp, aes(y = means, x = temperature)) +
  geom_ribbon(data = p.AF.w.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.AF.w.u.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  geom_line(data = p.AF.w.u.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "blue4") +
  
  geom_line(data = p.AF.m.u.l.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) +
  geom_ribbon(data = p.AF.m.u.l.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.AF.m.u.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.u.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.AF.m.u.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.u.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$max)), max(round(af_sig1$max)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 8) +
  ylab("") + 
  xlab("") +
  ggtitle("Direct Fire")

# IED:
plot.AF.m.i.l.qp <- ggplot() + 
  
  geom_line(data = p.AF.w.i.q.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  geom_line(data = p.AF.w.i.c.qp, aes(y = means, x = temperature)) + 
  geom_ribbon(data = p.AF.w.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "firebrick") +
  
  geom_line(data = p.AF.m.i.q.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.i.q.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  geom_line(data = p.AF.m.i.c.qp, aes(y = means, x = temperature), color = "gray20", alpha = 0.5) + 
  geom_ribbon(data = p.AF.m.i.c.qp, aes(x = temperature, y = means, ymin = cil, ymax = ciu), alpha = .15, fill = "gray20") +
  
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = seq(min(round(af_sig1$max)), max(round(af_sig1$max)), by = 6)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=12),
    axis.text.x = element_text(angle = 45)) + 
  ylim(-0.5, 8) +
  ylab("") + 
  xlab("") +
  ggtitle("IED Attacks")

# Combine Plots:

IQ.plots <- ggarrange(plot.IQ.m.u.l.qp + ylab("Iraq"), plot.IQ.m.c.l.qp, plot.IQ.m.i.l.qp, ncol = 3)
AF.plots <- ggarrange(plot.AF.m.u.l.qp + ylab("Afghanistan"), plot.AF.m.i.l.qp)

combined <- ggarrange(IQ.plots, AF.plots, nrow = 2)
combined <- annotate_figure(combined,
                            left = text_grob("Expected Attack Count", color = "black", rot = 90, size = 18),
                            bottom = text_grob("Temperature (°F)", color = "black", size = 18))

# Then, combine all figures:

# PRODUCES FIGURE 2 IN PAPER:
ggarrange(p_combined, combined)

###
### Baghdad Surveys
###

BG.m.l.b <- felm(vs_mnf ~ max +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + block + date_month | 0 | block, data = bg_sur)
BG.m.q.b <- felm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + block + date_month | 0 | block, data = bg_sur)
BG.m.c.b <- felm(vs_mnf ~ max +  max2 +  max3 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + block + date_month | 0 | block, data = bg_sur)
BG.m.l.m <- felm(vs_mnf ~ max +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + mahala + date_month | 0 | mahala, data = bg_sur)
BG.m.q.m <- felm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + mahala + date_month | 0 | mahala, data = bg_sur)
BG.m.c.m <- felm(vs_mnf ~ max +  max2 +  max3 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + daylight +
                   p3_electr_binary + 
                   temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 |
                   education + mahala + date_month | 0 | mahala, data = bg_sur)

stargazer(BG.m.q.b, BG.m.q.m)

BG.m.l.b.lg <- bayesglm(vs_mnf ~ max +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.q.b.lg <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.c.b.lg <- bayesglm(vs_mnf ~ max +  max2 +  max3 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(block) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))

BG.m.l.m.lg <- bayesglm(vs_mnf ~ max +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.q.m.lg <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))
BG.m.c.m.lg <- bayesglm(vs_mnf ~ max +  max2 +  max3 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family = binomial(link = "logit"))

#Next, repeat the analysis, using various income subsets:
bg_sur.s1.l <- bg_sur[bg_sur$income.w < 65000, ]
bg_sur.s2.l <- bg_sur[bg_sur$income.w < 70000, ]
bg_sur.s3.l <- bg_sur[bg_sur$income.w < 75000, ]
bg_sur.s4.l <- bg_sur[bg_sur$income.w < 80000, ]
bg_sur.s5.l <- bg_sur[bg_sur$income.w < 85000, ]
bg_sur.s1.u <- bg_sur[bg_sur$income.w > 65000, ]
bg_sur.s2.u <- bg_sur[bg_sur$income.w > 70000, ]
bg_sur.s3.u <- bg_sur[bg_sur$income.w > 75000, ]
bg_sur.s4.u <- bg_sur[bg_sur$income.w > 80000, ]
bg_sur.s5.u <- bg_sur[bg_sur$income.w > 85000, ]

BG.m.q.m.lg.s1.l <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s1.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s2.l <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s2.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s3.l <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s3.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s4.l <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s4.l, family = binomial(link = "logit"))
BG.m.q.m.lg.s5.l <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s5.l, family = binomial(link = "logit"))

BG.m.q.m.lg.s1.u <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s1.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s2.u <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s2.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s3.u <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s3.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s4.u <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s4.u, family = binomial(link = "logit"))
BG.m.q.m.lg.s5.u <- bayesglm(vs_mnf ~ max +  max2 +  age + work_wkhr + income.w + hhld_size + m_sig1_7 + prcp + wdsp + dewp + visib + mxspd + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur.s5.u, family = binomial(link = "logit"))

plot_input_function <- function(x, y){
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(x), Sigma = vcov(x))
  paste <- model.matrix(x)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps <- matrix(NA, length(round(min(bg_sur$max)):round(max(bg_sur$max))), 100)
  for(i in 1:length(round(min(bg_sur$max)):round(max(bg_sur$max)))){
    if(y %in% "linear"){
      paste[, 2] <- i + (round(min(bg_sur$max))-1)
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta) / (exp(x.beta) + 1)
      reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(y %in% "quadratic"){
      paste[, 2] <- i + (round(min(bg_sur$max))-1)
      paste[, 3] <- i + (round(min(bg_sur$max))-1)^2
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta) / (exp(x.beta) + 1)
      reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
    } else if(y %in% "cubic"){
      paste[, 2] <- i + (round(min(bg_sur$max))-1)
      paste[, 3] <- i + (round(min(bg_sur$max))-1)^2
      paste[, 4] <- i + (round(min(bg_sur$max))-1)^3
      x.beta <- paste %*% t(beta.sim)
      yhat <- exp(x.beta) / (exp(x.beta) + 1)
      reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
    }
  }
  return(reps)
}

# Next, given that the quadratic specification is the most significant, use that to 
# compare results across different income levels:

primary <- plot_input_function(BG.m.l.m.lg, "linear")
means.primary_l <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_l <- means.primary_l + qnorm(0.975)*se.primary
ciu.primary_l <- means.primary_l - qnorm(0.975)*se.primary
primary.dataset <- cbind(as.data.frame(means.primary_l), as.data.frame(cil.primary_l), as.data.frame(ciu.primary_l))

primary <- plot_input_function(BG.m.q.m.lg, "quadratic")
means.primary_q_m <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_q_m <- means.primary_q_m + qnorm(0.975)*se.primary
ciu.primary_q_m <- means.primary_q_m - qnorm(0.975)*se.primary
primary.dataset <- cbind(primary.dataset, as.data.frame(means.primary_q_m), as.data.frame(cil.primary_q_m), as.data.frame(ciu.primary_q_m))

primary <- plot_input_function(BG.m.q.b.lg, "quadratic")
means.primary_q_b <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_q_b <- means.primary_q_b + qnorm(0.975)*se.primary
ciu.primary_q_b <- means.primary_q_b - qnorm(0.975)*se.primary
primary.dataset <- cbind(primary.dataset, as.data.frame(means.primary_q_b), as.data.frame(cil.primary_q_b), as.data.frame(ciu.primary_q_b))

primary <- plot_input_function(BG.m.c.m.lg, "cubic")
means.primary_c <- apply(primary, 1 , mean) 
se.primary <- apply(primary, 1 , sd)
cil.primary_c <- means.primary_c + qnorm(0.975)*se.primary
ciu.primary_c <- means.primary_c - qnorm(0.975)*se.primary
primary.dataset <- cbind(primary.dataset, as.data.frame(means.primary_c), as.data.frame(cil.primary_c), as.data.frame(ciu.primary_c))

s1.l <- plot_input_function(BG.m.q.m.lg.s1.l, "quadratic")
means.s1.l <- apply(s1.l, 1 , mean) 
se.s1.l <- apply(s1.l, 1 , sd)
cil.s1.l <- means.s1.l + qnorm(0.975)*se.s1.l
ciu.s1.l <- means.s1.l - qnorm(0.975)*se.s1.l

s2.l <- plot_input_function(BG.m.q.m.lg.s2.l, "quadratic")
means.s2.l <- apply(s2.l, 1 , mean) 
se.s2.l <- apply(s2.l, 1 , sd)
cil.s2.l <- means.s2.l + qnorm(0.975)*se.s2.l
ciu.s2.l <- means.s2.l - qnorm(0.975)*se.s2.l

s3.l <- plot_input_function(BG.m.q.m.lg.s3.l, "quadratic")
means.s3.l <- apply(s3.l, 1 , mean) 
se.s3.l <- apply(s3.l, 1 , sd)
cil.s3.l <- means.s3.l + qnorm(0.975)*se.s3.l
ciu.s3.l <- means.s3.l - qnorm(0.975)*se.s3.l

s4.l <- plot_input_function(BG.m.q.m.lg.s4.l, "quadratic")
means.s4.l <- apply(s4.l, 1 , mean) 
se.s4.l <- apply(s4.l, 1 , sd)
cil.s4.l <- means.s4.l + qnorm(0.975)*se.s4.l
ciu.s4.l <- means.s4.l - qnorm(0.975)*se.s4.l

s5.l <- plot_input_function(BG.m.q.m.lg.s5.l, "quadratic")
means.s5.l <- apply(s5.l, 1 , mean) 
se.s5.l <- apply(s5.l, 1 , sd)
cil.s5.l <- means.s5.l + qnorm(0.975)*se.s5.l
ciu.s5.l <- means.s5.l - qnorm(0.975)*se.s5.l

s1.u <- plot_input_function(BG.m.q.m.lg.s1.u, "quadratic")
means.s1.u <- apply(s1.u, 1 , mean) 
se.s1.u <- apply(s1.u, 1 , sd)
cil.s1.u <- means.s1.u + qnorm(0.975)*se.s1.u
ciu.s1.u <- means.s1.u - qnorm(0.975)*se.s1.u

s2.u <- plot_input_function(BG.m.q.m.lg.s2.u, "quadratic")
means.s2.u <- apply(s2.u, 1 , mean) 
se.s2.u <- apply(s2.u, 1 , sd)
cil.s2.u <- means.s2.u + qnorm(0.975)*se.s2.u
ciu.s2.u <- means.s2.u - qnorm(0.975)*se.s2.u

s3.u <- plot_input_function(BG.m.q.m.lg.s3.u, "quadratic")
means.s3.u <- apply(s3.u, 1 , mean) 
se.s3.u <- apply(s3.u, 1 , sd)
cil.s3.u <- means.s3.u + qnorm(0.975)*se.s3.u
ciu.s3.u <- means.s3.u - qnorm(0.975)*se.s3.u

s4.u <- plot_input_function(BG.m.q.m.lg.s4.u, "quadratic")
means.s4.u <- apply(s4.u, 1 , mean) 
se.s4.u <- apply(s4.u, 1 , sd)
cil.s4.u <- means.s4.u + qnorm(0.975)*se.s4.u
ciu.s4.u <- means.s4.u - qnorm(0.975)*se.s4.u

s5.u <- plot_input_function(BG.m.q.m.lg.s5.u, "quadratic")
means.s5.u <- apply(s5.u, 1 , mean) 
se.s5.u <- apply(s5.u, 1 , sd)
cil.s5.u <- means.s5.u + qnorm(0.975)*se.s5.u
ciu.s5.u <- means.s5.u - qnorm(0.975)*se.s5.u

primary.dataset <- cbind(primary.dataset, 
                         as.data.frame(means.s1.l), as.data.frame(cil.s1.l), as.data.frame(ciu.s1.l),
                         as.data.frame(means.s2.l), as.data.frame(cil.s2.l), as.data.frame(ciu.s2.l),
                         as.data.frame(means.s3.l), as.data.frame(cil.s3.l), as.data.frame(ciu.s3.l),
                         as.data.frame(means.s4.l), as.data.frame(cil.s4.l), as.data.frame(ciu.s4.l),
                         as.data.frame(means.s5.l), as.data.frame(cil.s5.l), as.data.frame(ciu.s5.l),
                         as.data.frame(means.s1.u), as.data.frame(cil.s1.u), as.data.frame(ciu.s1.u), 
                         as.data.frame(means.s2.u), as.data.frame(cil.s2.u), as.data.frame(ciu.s2.u),
                         as.data.frame(means.s3.u), as.data.frame(cil.s3.u), as.data.frame(ciu.s3.u),
                         as.data.frame(means.s4.u), as.data.frame(cil.s4.u), as.data.frame(ciu.s4.u),
                         as.data.frame(means.s5.u), as.data.frame(cil.s5.u), as.data.frame(ciu.s5.u))

# Next, bound all negative predicted probabilities at 0 and 1:
primary.dataset <- primary.dataset %>% mutate(across(everything(), ~replace(., . < 0 , 0)))
primary.dataset <- primary.dataset %>% mutate(across(everything(), ~replace(., . > 1 , 1)))
# For plotting purposes, multiply all values by 100:
primary.dataset <- primary.dataset * 100

primary.dataset$temperature <-  round(min(bg_sur$max)):round(max(bg_sur$max))

distinct_survey_days <- distinct(bg_sur, date, .keep_all = TRUE)

p1 <- ggplot(data=primary.dataset) +

  geom_line(aes(x=temperature, y=means.primary_q_b), colour='blue4') +
  geom_ribbon(aes(x=temperature, ymin = cil.primary_q_b, ymax = ciu.primary_q_b), alpha = .2, fill = "blue4") +
  
  geom_line(aes(x=temperature, y=means.primary_q_m), colour='gray', alpha = 1) +
  geom_ribbon(aes(x=temperature, ymin = cil.primary_q_m, ymax = ciu.primary_q_m), alpha = .2, fill = "gray") +

theme_bw() +
  ggtitle("Full Iraqi Male Sample") +
  theme(text=element_text(size=14), 
        legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), 
        plot.background=element_blank(), 
        axis.title.x=element_blank(),
        axis.title.y=element_blank()) +
  ylim(40,100) +
  geom_vline(xintercept = round(max(bg_sur$max)), linetype=4) +
  xlim(round(min(bg_sur$max)), round(max(bg_sur$max))) + 
  scale_y_continuous(breaks = seq(40,100, by = 20), labels = scales::percent(seq(0.4,1, by = 0.20))) + 
  scale_x_continuous(breaks = seq(round(min(bg_sur$max)), round(max(bg_sur$max)), by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=14),
    axis.text.x = element_text(angle = 45))

p2 <- ggplot(data=primary.dataset) +
  geom_line(aes(x=temperature, y=means.s1.l, colour="< 65,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s2.l, colour="< 70,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s3.l, colour="< 75,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s4.l, colour="< 80,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s5.l, colour="< 85,000 IQD"), linetype="solid", alpha = 0.5) +
  
  geom_ribbon(aes(x=temperature, ymin = cil.s1.l, ymax = ciu.s1.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s2.l, ymax = ciu.s2.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s3.l, ymax = ciu.s3.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s4.l, ymax = ciu.s4.l), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s5.l, ymax = ciu.s5.l), alpha = .1, fill = "blue4") +
  
  ylim(32.5, 100) +
  
  theme_bw() +
  ggtitle("Highest Incomes Excluded") +
  scale_colour_manual("", 
                      breaks = c("< 65,000 IQD", "< 70,000 IQD", "< 75,000 IQD", "< 80,000 IQD", "< 85,000 IQD"),
                      values = c("firebrick", "green4", "blue4", "yellow4", "orange4")) +
  theme(text=element_text(size=14), 
        legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), 
        plot.background=element_blank(), 
        axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_vline(xintercept = round(max(bg_sur$max)), linetype=4) +
  xlim(round(min(bg_sur$max)), round(max(bg_sur$max))) + 
  scale_y_continuous(breaks = seq(0,100, by = 20), labels = scales::percent(seq(0,1, by = 0.20))) + 
  scale_x_continuous(breaks = seq(round(min(bg_sur$max)), round(max(bg_sur$max)), by = 5)) + 
  theme(
    axis.ticks = element_blank(), 
    text = element_text(size=14),
    axis.text.x = element_text(angle = 45))

p3 <- ggplot(data=primary.dataset) +
  geom_histogram(data = distinct_survey_days, aes(x=max, y=(..count..)), bins = 100, size=.1, color="black", alpha=.4) +
  
  geom_line(aes(x=temperature, y=means.s1.u, colour="> 65,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s2.u, colour="> 70,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s3.u, colour="> 75,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s4.u, colour="> 80,000 IQD"), linetype="solid", alpha = 0.5) +
  geom_line(aes(x=temperature, y=means.s5.u, colour="> 85,000 IQD"), linetype="solid", alpha = 0.5) +
  
  geom_ribbon(aes(x=temperature, ymin = cil.s1.u, ymax = ciu.s1.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s2.u, ymax = ciu.s2.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s3.u, ymax = ciu.s3.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s4.u, ymax = ciu.s4.u), alpha = .1, fill = "blue4") +
  geom_ribbon(aes(x=temperature, ymin = cil.s5.u, ymax = ciu.s5.u), alpha = .1, fill = "blue4") +
  ylim(0, 100) +
  theme_bw() +
  ggtitle("Highest Incomes Only") +
  scale_colour_manual("", 
                      breaks = c("> 65,000 IQD", "> 70,000 IQD", "> 75,000 IQD", "> 80,000 IQD", "> 85,000 IQD"),
                      values = c("firebrick", "green4", "blue4", "yellow4", "orange4")) +
  theme(text=element_text(size=14), 
        legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), 
        plot.background=element_blank(), 
        axis.title.x=element_blank(),
        axis.title.y=element_blank())+
  geom_vline(xintercept = round(max(bg_sur$max)), linetype=4) +
  xlim(round(min(bg_sur$max)), round(max(bg_sur$max))) + 
  scale_y_continuous(breaks = seq(0,100, by = 20), labels = scales::percent(seq(0,1, by = 0.20))) + 
  scale_x_continuous(breaks = seq(round(min(bg_sur$max)), round(max(bg_sur$max)), by = 5)) + 
  theme(
    #axis.text.y = element_blank(), 
    axis.ticks = element_blank(), 
    text = element_text(size=14),
    axis.text.x = element_text(angle = 45))

# PRODUCES FIGURE 5 IN PAPER:
plot.attitudes <- grid.arrange(p1, arrangeGrob(p2, p3, ncol=1), 
                               ncol=2, widths=c(1,1.2), bottom = textGrob("Temperature (°F)", gp=gpar(fontsize=18)), 
                               left = textGrob("Predicted Probability", rot = 90, gp=gpar(fontsize=18)), 
                               right = "", 
                               top="")

###
### END.
###

###
### COMPARISONS:
###

# Next, compare the day-level temperature results to other factors associated 
# with differing levels of insurgent violence or support for violence:

# First, how does support for violence vary across those with differing
# perceptions of past conditions?

# Create an indicator of average perceptions across variables:
bg_sur$perceptions <- apply(bg_sur[, c("p3_fam_binary", "p3_govteff_binary", "p3_jobs_binary", 
                                       "p3_electr_binary", "p3_sec_binary", "p3_iping_binary", 
                                       "p3_baghdad_binary")], 1, mean, na.rm = T)

glm.sat.quad.p3_fam <- glm(vs_mnf ~ p3_fam_binary +
                             age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_govteff <- glm(vs_mnf ~ p3_govteff_binary + 
                                 age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_jobs <- glm(vs_mnf ~ p3_jobs_binary + 
                              age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_electr <- glm(vs_mnf ~ p3_electr_binary + 
                                age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_sec <- glm(vs_mnf ~ p3_sec_binary + 
                             age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_iping <- glm(vs_mnf ~ p3_iping_binary + 
                               age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_baghdad <- glm(vs_mnf ~ p3_baghdad_binary + 
                                 age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))
glm.sat.quad.p3_perceptions <- glm(vs_mnf ~ perceptions +
                                     age + work_wkhr + income.w + hhld_size + m_sig1_7 + as.factor(education) + as.factor(mahala) + as.factor(date_month), data = bg_sur, family=binomial(link="logit"))

# Export regression results:
stargazer(glm.sat.quad.p3_fam,
          glm.sat.quad.p3_govteff,
          glm.sat.quad.p3_jobs,
          glm.sat.quad.p3_electr,
          glm.sat.quad.p3_sec,
          glm.sat.quad.p3_iping,
          glm.sat.quad.p3_baghdad,
          glm.sat.quad.p3_perceptions,
          title = "",
          model.names = TRUE,
          omit.yes.no = c("Yes", "No"),
          star.cutoffs = c(.1, .05, .01),
          omit = "factor", single.row = F)

# Generate predicted probabilities:
plot_input_function_binary <- function(x){
  set.seed(1)
  beta.sim <- mvrnorm(100, mu = coef(x), Sigma = vcov(x))
  paste <- model.matrix(x)
  names1 <- colnames(paste)
  beta.sim <- as.data.frame(beta.sim)
  beta.sim <- beta.sim[, c(names1)]
  
  reps <- matrix(NA, length(0:1), 100)
  for(i in 1:length(0:1)){
    paste[, 2] <- i - 1 
    x.beta <- paste %*% t(beta.sim)
    yhat <- exp(x.beta) / (exp(x.beta) + 1)
    reps[i, ] <- apply(yhat, 2, mean, na.rm=T)
  }
  return(reps)
}

# For past perceptions:
p3_1.u <- plot_input_function_binary(glm.sat.quad.p3_fam)
means.p3_1.u <- apply(p3_1.u, 1 , mean)
se.p3_1.u <- apply(p3_1.u, 1 , sd)
cil.p3_1.u <- means.p3_1.u + qnorm(0.975)*se.p3_1.u
ciu.p3_1.u <- means.p3_1.u - qnorm(0.975)*se.p3_1.u

p3.dataset <- cbind(as.data.frame(means.p3_1.u), as.data.frame(cil.p3_1.u), as.data.frame(ciu.p3_1.u), as.data.frame(rep("X",2)))
colnames(p3.dataset)[4] <- "X"
p3.dataset$p3 <- 0:1

p3_2.u <- plot_input_function_binary(glm.sat.quad.p3_govteff)
means.p3_2.u <- apply(p3_2.u, 1 , mean)
se.p3_2.u <- apply(p3_2.u, 1 , sd)
cil.p3_2.u <- means.p3_2.u + qnorm(0.975)*se.p3_2.u
ciu.p3_2.u <- means.p3_2.u - qnorm(0.975)*se.p3_2.u

p3_3.u <- plot_input_function_binary(glm.sat.quad.p3_jobs)
means.p3_3.u <- apply(p3_3.u, 1 , mean)
se.p3_3.u <- apply(p3_3.u, 1 , sd)
cil.p3_3.u <- means.p3_3.u + qnorm(0.975)*se.p3_3.u
ciu.p3_3.u <- means.p3_3.u - qnorm(0.975)*se.p3_3.u

p3_4.u <- plot_input_function_binary(glm.sat.quad.p3_electr)
means.p3_4.u <- apply(p3_4.u, 1 , mean)
se.p3_4.u <- apply(p3_4.u, 1 , sd)
cil.p3_4.u <- means.p3_4.u + qnorm(0.975)*se.p3_4.u
ciu.p3_4.u <- means.p3_4.u - qnorm(0.975)*se.p3_4.u

p3_5.u <- plot_input_function_binary(glm.sat.quad.p3_sec)
means.p3_5.u <- apply(p3_5.u, 1 , mean)
se.p3_5.u <- apply(p3_5.u, 1 , sd)
cil.p3_5.u <- means.p3_5.u + qnorm(0.975)*se.p3_5.u
ciu.p3_5.u <- means.p3_5.u - qnorm(0.975)*se.p3_5.u

p3_6.u <- plot_input_function_binary(glm.sat.quad.p3_iping)
means.p3_6.u <- apply(p3_6.u, 1 , mean)
se.p3_6.u <- apply(p3_6.u, 1 , sd)
cil.p3_6.u <- means.p3_6.u + qnorm(0.975)*se.p3_6.u
ciu.p3_6.u <- means.p3_6.u - qnorm(0.975)*se.p3_6.u

p3_7.u <- plot_input_function_binary(glm.sat.quad.p3_baghdad)
means.p3_7.u <- apply(p3_7.u, 1 , mean)
se.p3_7.u <- apply(p3_7.u, 1 , sd)
cil.p3_7.u <- means.p3_7.u + qnorm(0.975)*se.p3_7.u
ciu.p3_7.u <- means.p3_7.u - qnorm(0.975)*se.p3_7.u

p3_8.u <- plot_input_function_binary(glm.sat.quad.p3_perceptions)
means.p3_8.u <- apply(p3_8.u, 1 , mean)
se.p3_8.u <- apply(p3_8.u, 1 , sd)
cil.p3_8.u <- means.p3_8.u + qnorm(0.975)*se.p3_8.u
ciu.p3_8.u <- means.p3_8.u - qnorm(0.975)*se.p3_8.u

p3.dataset <- cbind(p3.dataset,
                    as.data.frame(means.p3_2.u), as.data.frame(cil.p3_2.u), as.data.frame(ciu.p3_2.u),
                    as.data.frame(means.p3_3.u), as.data.frame(cil.p3_3.u), as.data.frame(ciu.p3_3.u),
                    as.data.frame(means.p3_4.u), as.data.frame(cil.p3_4.u), as.data.frame(ciu.p3_4.u),
                    as.data.frame(means.p3_5.u), as.data.frame(cil.p3_5.u), as.data.frame(ciu.p3_5.u),
                    as.data.frame(means.p3_6.u), as.data.frame(cil.p3_6.u), as.data.frame(ciu.p3_6.u), 
                    as.data.frame(means.p3_7.u), as.data.frame(cil.p3_7.u), as.data.frame(ciu.p3_7.u),
                    as.data.frame(means.p3_8.u), as.data.frame(cil.p3_8.u), as.data.frame(ciu.p3_8.u))

# PRODUCES FIGURE 20 IN PAPER:
p3 <- ggplot(data=p3.dataset) +
  geom_point(aes(x=p3, y=means.p3_1.u, colour="Family")) +
  geom_point(aes(x=p3, y=means.p3_2.u, colour="Government")) +
  geom_point(aes(x=p3, y=means.p3_3.u, colour="Jobs")) +
  geom_point(aes(x=p3, y=means.p3_4.u, colour="Electricity")) +
  geom_point(aes(x=p3, y=means.p3_5.u, colour="Security")) +
  geom_point(aes(x=p3, y=means.p3_6.u, colour="Iraqi Police")) +
  geom_point(aes(x=p3, y=means.p3_7.u, colour="Baghdad")) +
  geom_point(aes(x=p3, y=means.p3_8.u, colour="Average")) +
  
  geom_errorbar(aes(x=p3, ymin=cil.p3_1.u,ymax=ciu.p3_1.u,colour="Family",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_2.u,ymax=ciu.p3_2.u,colour="Government",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_3.u,ymax=ciu.p3_3.u,colour="Jobs",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_4.u,ymax=ciu.p3_4.u,colour="Electricity",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_5.u,ymax=ciu.p3_5.u,colour="Security",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_6.u,ymax=ciu.p3_6.u,colour="Iraqi Police",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_7.u,ymax=ciu.p3_7.u,colour="Baghdad",width=0.1)) +
  geom_errorbar(aes(x=p3, ymin=cil.p3_8.u,ymax=ciu.p3_8.u,colour="Average",width=0.1)) +
  
  theme_bw() +
  ggtitle("") +
  xlab("Negative Perceptions Indicator") +
  ylab("Predicted Probability") +
  theme(text=element_text(size=26),
        legend.title=element_blank(),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank()) +
  scale_x_continuous(breaks = c(0, 1))
p3

# To what extent to levels of violence in Baghdad and Basrah observed during the study period differ as a
# matter of the day of week?
# Create day of week indicator:
iq_sig1$day_of_week <- as.factor(wday(iq_sig1$date))
# Set Fridays, the day of prayer, to the omitted comparison:
iq_sig1$day_of_week <- relevel(iq_sig1$day_of_week, ref=6)
# Create a measure of all insurgent violence:
iq_sig1$all_viol <- iq_sig1$unconstrained + iq_sig1$constrained + iq_sig1$ied + iq_sig1$idf

nb.all_viol.dow <- glm.nb(all_viol ~ as.factor(day_of_week) + temp + temp2 + prcp + wdsp + dewp + visib + mxspd + daylight + hoursofpower +
                            temp_1 + temp_2 + temp_3 + temp_4 + temp_5 + temp_6 + temp_7 +
                            unconstrained_1 + unconstrained_2 + unconstrained_3 + unconstrained_4 + unconstrained_5 + unconstrained_6 + unconstrained_7 +
                            constrained_1 + constrained_2 + constrained_3 + constrained_4 + constrained_5 + constrained_6 + constrained_7 +
                            ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 +
                            as.factor(governorate) + as.factor(week), data = iq_sig1)

# Plot the expected counts of insurgent violence by day of week:
set.seed(1)
beta.sim <- mvrnorm(1000, mu = coef(nb.all_viol.dow), Sigma = vcov(nb.all_viol.dow))
paste <- model.matrix(nb.all_viol.dow)
names1 <- colnames(paste)
beta.sim <- as.data.frame(beta.sim)
beta.sim <- beta.sim[, c(names1)]

reps.dow <- matrix(NA, length(1:7), 1000)
for(i in 1:length(1:7)){
  if(i == 1){
    paste[, 2:7] <- 0
    paste[, 2] <- 1
  } else if(i == 2){
    paste[, 2:7] <- 0
    paste[, 3] <- 1
  } else if(i == 3){
    paste[, 2:7] <- 0
    paste[, 4] <- 1
  } else if(i == 4){
    paste[, 2:7] <- 0
    paste[, 5] <- 1
  } else if(i == 5){
    paste[, 2:7] <- 0
    paste[, 6] <- 1
  } else if(i == 6){
    paste[, 2:7] <- 0 # to capture base.
  } else if(i == 7){
    paste[, 2:7] <- 0
    paste[, 7] <- 1
  }
  x.beta <- paste %*% t(beta.sim)
  yhat <- exp(x.beta)
  reps.dow[i, ] <- apply(yhat, 2, mean, na.rm=T)
}

means.dow <- apply(reps.dow, 1 , mean) 
se.dow <- apply(reps.dow, 1 , sd)
cil.dow <- means.dow + qnorm(0.975)*se.dow
ciu.dow <- means.dow - qnorm(0.975)*se.dow

# PRODUCES FIGURE 19 IN PAPER:
plot.dow <- qplot(1:7, means.dow, geom=c("point"), ylim = c(15, 22.5), ylab = "", xlab = "", main = "All Insurgent Violence", colour=I("red")) + 
  geom_errorbar(ymin = cil.dow, ymax = ciu.dow, alpha = .2) +
  theme_bw() +
  theme(text = element_text(size=26), legend.title=element_blank(), 
        panel.background=element_blank(),
        panel.grid.minor=element_blank(), plot.background=element_blank()) +
  theme(legend.key = element_blank()) +
  scale_x_continuous(breaks = c(1, 2, 3, 4, 5, 6, 7), labels = c("Sun", "Mon", "Tues", "Wed", "Thur", "Fri", "Sat"))
plot.dow

###
### END.
###
